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Abstract 

Modeling relations between individuals is a classical question in social sciences 
and clustering individuals according to the observed patterns of interactions allows 
to uncover a latent structure in the data. Stochastic block model (SBM) is a popu¬ 
lar approach for grouping the individuals with respect to their social comportment. 
When several relationships of various types can occur jointly between the individu¬ 
als, the data are represented by multiplex networks where more than one edge can 
exist between the nodes. In this paper, we extend the SBM to multiplex networks 
in order to obtain a clustering based on more than one kind of relationship. We pro¬ 
pose to estimate the parameters -such as the marginal probabilities of assignment 
to groups (blocks) and the matrix of probabilities of connections between groups- 
* pierre.barbillon@agroparistech.fr 
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through a variational Expectation-Maximization procedure. Consistency of the es¬ 
timates as well as statistical properties of the model are obtained. The number of 
groups is chosen thanks to the Integrated Completed Likelihood criteria, a penalized 
likelihood criterion. Multiplex Stochastic Block Model arises in many situations but 
our applied example is motivated by a network of French cancer researchers. The 
two possible links (edges) between researchers are a direct connection or a connec¬ 
tion through their labs. Our results show strong interactions between these two 
kinds of connections and the groups that are obtained are discussed to emphasize 
the common features of researchers grouped together. 

Keywords — Bivariate Stochastic Block Model, Multilevel / Multiplex networks, 
Social network. 

1 Introduction 

Network analysis has emerged as a key technique for understanding and for investigating 
social interactions through the properties of relations between and within units. From a 
statistical point of view, a network is a realization of a random graph formed by a set of 
nodes V representing the units (e.g. individuals, actors, companies) and a set of edges 
E representing relationships between pairs of nodes. 

The system in which the same nodes belong to multiple networks is typically referred 
to as a multiplex network or multigraph (see Wasserman (1994) for example). In recent 
literature, there has been an upsurge of interest in multiplex networks (see for example 
Cozzo et al. (2012); Loe and Jeldtoft Jensen (2014); Rank et al. (2010); Szell et al. (2010); 
Mucha et al. (2010); Maggioni et al. (2013); Brummitt et al. (2012); Saumell-Mendiola 
et al. (2012); Bianconi (2013); Nicosia et al. (2013)). In these multiplex networks, differ¬ 
ent kinds of links (or connections) are possible for each pair of nodes. This induced link 
multiplexity is a fundamental aspect of social relations (Snijders and Baerveldt, 2003) 
since these multiple links are frequently interdependent: links in one network may have 
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an influence on the formation or dissolution of links in other networks. 

The simultaneous analysis of several networks also arises when one is interested in the 
social comportment of individuals belonging to organized entities (such as companies, 
laboratories, political groups, etc.), with some individuals possibly belonging to the same 
institution. While the actors will exchange resources (such as advice for instance) at the 
individual level, their respective organizations of affiliation will also share resources at 
the institutional level (financial resources for instance). Each level (individuals and or¬ 
ganizations) constitutes a system of exchange of different resources that has its own 
logic and could be studied separately. However, studying the two networks jointly (and 
hence embedding the individuals in the multilevel relational and organizational struc¬ 
tures constituting the inter-organizational context of their actions) would allow us to 
identify the individuals that benefit from relatively easy access to the resources circu¬ 
lating in each level, which is of much more interest. In other words, studying the two 
levels jointly could help us understand how an individual can benefit from the position 
of its organization in the institutional network. 

In this paper, we are interested in studying the advice relations between researchers 
and the exchanges of resources between laboratories. We adopt the following individual- 
oriented strategy (this point is discussed in the paper): the institutional network is used 
to define a new network on the individual level i.e. the set of nodes consists in the set 
of individuals and for a pair of individuals, two kinds of links are possible: a direct con¬ 
nection given by the individual network and a connection through their organizations 
given by the organizational network. As a consequence, the individual and institutional 
levels are fused into a multiplex. 

We then develop a statistical model able to detect in multiplex substantial non¬ 
trivial topological features, with patterns of connection between their elements that are 
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not purely regular. Several models such as scale-free networks and small-world networks 
have been proposed to describe and understand the heterogeneity observed in networks. 
These models allow to derive properties of the network at the macro-scale and to un¬ 
derstand the outcomes of interactions. To explore heterogeneity at others scales (such 
as micro or meso-scale) in social networks, specific models such as the stochastic block 
models (SBM) (Snijders and Nowicki (1997)) have been developed for uniplex networks. 
In this paper, we propose an original extension of the SBMs to the multiplex case. Our 
model is efficient to model not only the main effects (that correspond to a classical uni¬ 
plex) but also the pairwise interactions between the nodes. We estimate the parameters 
of the multiplex SBM using an extension of the variational EM algorithm. Consistency 
of the estimation of the parameters is proved. As for uniplex SBM, a key issue is to 
choose the number of blocks. We use a penalized likelihood criterion, namely Integrated 
Completed Likelihood (ICL). The inference procedure is performed on the cancer re¬ 
searchers / laboratories dataset. 

The paper is organized as follows. The extension of SBM to multiplex network is pre¬ 
sented in Section 2, the proofs of model identifiability and the consistency of variational 
EM procedure are postponed in Appendices A and B. In Section 3, we describe Lazega 
et al. (2008) ’s dataset, apply the new modeling and discuss the results. Eventually, the 
contribution of multiplex SBM to the analysis of multiplex networks is highlighted in 
Section 4. 

2 Multiplex stochastic block model 

The main objective is to cluster the individuals (or nodes) into blocks sharing connection 
properties with the other individuals of the multiplex-network. Stochastic block models 
(Nowicki and Snijders, 2001) for random graphs have emerged as a natural tool to 
perform such a clustering based on uniplex networks (directed or not, valued or not). 
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In the following, we propose an extension of the Stochastic Block Model (SBM) to 
multiplex networks. The SBM for multiplex networks is derived from a multiplex Erdos- 
Renyi model which is described in subsection 2.1. The SBM for multiplex networks is 
derived in subsection 2 . 2 . 


2.1 Erdos-Renyi model for multiplex networks 

Let X 1 ,...,X A be K directed graphs relying on the same set of nodes E = {1,... ,n}. 
We assume that V(i,j),i / j, V7c £ {1,..., K}, £ {0,1} and Xu ^ 0. We define a 

joint distribution on X 1:A = (X 1 ,..., X A ) as: V(i, j) £ {1,..., n} 2 , i ^ j, \/w £ {0, 1} K , 


P (X}j K = w) = ttM where ^ w) = 1, 

ioe{o,i}^ 


( 1 ) 


and (Xjj K )i t] are mutually independent. 

The maximum likelihood estimate of the parameter of interest 7 r 
is, for all w £ { 0 , 1} K \ 


'K'W 


1 

n(n — 1 ) 


hxj: K =w}- 




This model is quite simple since any relation between two individuals (a relation being 
a collection of edges) does not depend on the relations between the other individuals. 
However, the different kind of relations between two individuals (edges) are not assumed 
to be independent. 


Remark 1. This model is clearly an extension of the Erdos-Renyi model since the 
marginal distribution of X(( (for any k = 1... K) is Bernoulli with density : 

«■ (*« = 4)=( X * < ” ,) ) ( E ' < " ) 

\™e{o,i}|«j fe =i / \M!e{o,i}|«) fe =o 

Moreover, any conditional distribution of Xf- given (Xk )/ £l 5 (where S\ k is a subset of 
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{1 not containing k) is also univariate Bernoulli. For instance, if K = 2 the 

conditional distribution of Xj- given Xf- is 


nx* = 4 x ?i = 4) = 




^’4) 


(!— 


ij i ■‘-ij 




(l,a:?.) , (0,x 2 .) 

7P 7P V' 


( 1 , 1 ?.) , (0,a;?.) 

7P W + 7P v' 


Moreover, the components of the bivariate Bernoulli random vector (Xh, X'F) are inde¬ 
pendent if and only if ir^ 00 ')= 7r( 10 )7r( 01 ). 


Introduction of explanatory variables. Naturally, the Erdos-Renyi model for mul¬ 
tiplex networks can be extended to take into account explanatory variables. Let y ^ 
denote the covariates characterizing the couple of nodes the model is defined by 

the probabilities: 


nxij K = v) 

P(X^ = (0,...,0)) 


exp( M (") + (/?(")) T yij ) 

1 +E„^(o,...,o) ex p(M (, ')+(/3 ( ’' ) ) T yo) 

_l_ 

1 +E^(o,...,o) ex p(/ iW +(^ w ) T yo) 


\/w / ( 0 ,..., 0 ), 


(2) 


where x T denotes the transposed vector of x. 

Remark 2. Note that in the multiplex Erdos-Renyi model, the modeling is actor based, 
which means that the individuals are the same for all the networks X 1 ,... ,X K and we 
model conjointly all the connections. As a consequence, the covariates y t j only depend 
on the couple ( i,j ) and are not linked to the network under consideration. 

Since the multiplex Erdos-Renyi model belongs to exponential models, the general¬ 
ized linear model theory applies when we introduce the covariates as in model ( 2 ) and 
the estimates are obtained using standard optimization strategies. 


2.2 Stochastic block model for multiplex networks 

When the goal is to cluster individuals according to their social comportment, we can 
derive a Stochastic Block Model version of the multiplex Erdos-Renyi model. Let Q be 
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the number of blocks and Z % the latent variable such that Z, = q if the individual i 
belongs to block q (note that an individual can only belong to one block in this version). 


The multiplex version of SBM is written as follows: V(?',j) G {l,...,n} 2 ,i / j, 
Vw <E {0,1} A , V(g, 0 G {1, • • •, Q} 2 , 


F(XljK = w\Z i = q,Z j =l ) = 

P(Z 2 ; = g) = a q . 

Such a model includes (2 K — 1)Q 2 + (Q — 1) parameters. Introducing the following 
notations: 


OL (cKl, . . . , CKq), TT (lTqi ^ (q, 7t) , 

the likelihood function is written as: 


£{X 1:K ]6) = f p(X 1:A '|Z;7r)p(Z;a)dZ, 

Jze{i,...Q} n 

Ze{l,...Q} n i,j,v£j *=1 

where the latent variable (the block affectations) are integrated out. The identifiability 
of the model can be proved (see Appendix A, theorem A.l) and the maximum likelihood 
estimators are consistent (theorem A.2). 


Remark 3. Note that, as before, covariates on the couple ( i,j ) can be introduced in the 
model: Vw ^ (0 ,..., 0), 


nXjf K = w\Zi = q, Zj = l) 
P(X^ = (0,...,0)| Zi = q,Zj=l) 


eX P (^i" ) + (/ 3 qr ) ) T yo) 
i+Ev^O.-.O) eX p(^I ) + (/ 3 9r ) ) Ty o) ’ 
_ 1 _ 

TEv/tO. 0 ) ex p(MqI ) + ( P ) T y a ) 


(5) 
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However, the number of parameters in this new model will increase drastically, leading 
to estimation issues. 

Maximum likelihood and model selection As soon as n or Q are large, the ob¬ 
served likelihood (4) is not tractable (due to the sum on Z 6 {1,... Q} n ) and its max¬ 
imization is a challenging task. Several approaches have been developed in the litera¬ 
ture (for a review, see Matias and Robin, 2014), both in the frequentist and Bayesian 
frameworks, starting from Snijders and Nowicki (1997) and Nowicki and Snijders (2001). 
However, when the latent data space is really large, these techniques can be burdensome. 
Some other strategies have been proposed, such as Bickel and Chen (2009) which relying 
on a profile-likelihood optimization or the moment estimation proposed by Ambroise 
and Matias (2012), to name but a few . 

The variational EM in the context of SBM proposed by Daudin et al. (2008) is a 
flexible tool to tackle the computational challenge in many types of graphs. Simulation 
studies showed its practical efficiency (Mariadassou et al., 2010). Moreover, its theoret¬ 
ical convergence towards the maximum likelihood estimates has been studied by Celisse 
et al. (2012) for binary graphs. In this paper, we adapt the variational EM to multiplex 
networks. The algorithm is described in Appendix B and its convergence towards the 
true parameter is proved (Theorem B.l). 

The selection of the most adequate number of blocks Q is performed using a mod¬ 
ification of the ICL criterion (as in Daudin et al., 2008; Mariadassou et al., 2010). Let 
A {q denote the model defined (5) with Q blocks. 

ICL(Mq) = maxlogp(X 1:X ,Z;60 - \ {Q 2 (2 k - l)\og(Kn(n- 1)) + (Q — l)logn} . 
6 Z 

where Z are the predictions of the assignments Z obtained as a sub-product of the 
variational EM algorithm (see Appendix B). As in the BIC criteria, the log refers to 
the number of data. Thus, the n nodes are used to estimate the Q — 1 probabilities 
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cki, ..., olq- i- The I\n{n — 1) edges are used to estimate 7r. No theoretical results exist 
for the ICL properties, but this criteria has proved its efficiency in practice. 

Remark 4. If we consider covariates as in equation (5), then the ICL is adapted : 
ICL(Mq) = max logp(X 1:A , Z; a, /j, (3) - ^ {Pq \og(Kn(n - 1)) + (Q - 1) logn} , 

aL,n,f3 2 

where Pq is the dimension of (/r, ( 3 ). 

3 Analysis for the laboratory-researcher data 

3.1 The data 

French scandals during the 1990s involving the voluntary sector around the cancer re¬ 
search dried up large donations that funded research laboratories. In the 2000s, the 
cancer research became politicized, with the launch of the Cancer Plan and the creation 
of a dedicated institution. The aim of this public agency is to coordinate the cancer 
research and to promote collaborations about top researchers. In this context, Lazega 
et al. (2008) studied the relations of advice between French cancer researchers identified 
as “Elite” conjointly with the relations of their respective laboratories. 

At the inter-individual level, the actors (researchers) were submitted a list of cancer 
researchers and asked in interviews whom they sought advice from. The advice were of 
five types, namely advice to deal with choices about the direction of projects, advice to 
find institutional support, advice to handle financial resources, advice for recruitment, 
and finally advice about manuscripts before submitting them to journals. The five advice 
networks are too sparse to be studied separately. That is why they were aggregated. 
Therefore two researchers are considered as linked if at least one kind of relationship 
exists. Obviously the links are directed. 

At the laboratory level (concerning only laboratories with “Elite” researchers), the lab- 
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oratory directors were asked to specify what type of resources they exchanged with the 
other laboratories on the list. The examined resources were the recruitment of post-docs 
and researchers, the development of programs of joint research, joint responses to ten¬ 
der offers, sharing of technical equipment, sharing of experimental material, mobility of 
administrative personnel, and invitations to conferences and seminars. Once again, to 
avoid over-sparsity, the various networks were aggregated and two laboratories are said 
to be linked if there is at least one link between them. From this network on labs, we 
can derive indirect links between the researchers, i.e. two researchers are connected if 
their laboratories exchange resources. We finally have two adjacency matrices on the 
same set of nodes (researchers). 

This corresponds to transforming the multilevel network (individual/ organization) 
into a multiplex network (several kinds of relations among individuals). This is rea¬ 
sonable since the majority of laboratories contains a unique “Elite” researcher. Thus, 
there is no big difference in the number of nodes between the institutional and individual 
levels. 

In addition, auxiliary covariates are available to describe the researchers: their age, 
their specialty, two publication performance scores based on two periods of five years, 
their status (director of the lab or not). Auxiliary covariates are also available the lab¬ 
oratories: their size (number of researchers) and their location. 

Complete data for 95 researchers identified as the “Elite” of French cancer research 
working in 76 laboratories are available. 


3.2 Statistical inference through SBM for multiplex 

To estimate the parameters of the multiplex SBM model on this dataset, we use a 
modified version of the variational EM algorithm (Daudin et ah, 2008) described in the 
Appendix Section B. The optimization of the ICL criterion derived from likelihood (4) 
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leads to four blocks (indistinctly denominated clusters or groups). For the sake of clarity, 
we index by R and L (rather than X 1 , X 2 ) the two adjacency matrices (respectively the 
direct and indirect ones) 

In Figures 1 and 2 we plot the marginal and conditional probabilities of the connec¬ 
tions of researchers (respectively labs) between and within blocks. 

Note that the study of the estimated marginal distributions allows us to have results 
on the researchers without considering the laboratories. This gives a clear interpreta¬ 
tion of the importance of the lab for the researcher network structure. The obtained 
blocks are described in Table 1: (la) gives the sizes of the four blocks; (lb), (lc) and 
(Id) describe the blocks with respect to the covariates “location”, “director of not” 
and “specialty”. The estimations by the variational EM procedure were conducted by 
wmixnet (Leger, 2014) with our specific implementation of the bivariate Bernoulli model. 

We now discuss the results. Multiplex SBM reveals interesting structural features of 
the multiplex network. More precisely, collaboration takes place in a clustered manner 
for both researchers and laboratories; collaborating laboratories tend to have affiliated 
researchers seeking advice from one another. Indeed, Figure 1 shows that the existence 
of a connection (exchange of resources) between labs clearly increases the probability of 
connection (sharing advice) between researchers. The reinforcement of this probability 
of connection is clearly outstanding in block 2. In this block, the researcher connections 
are quite unlikely within the block or with other blocks. However, conditionally to the 
existence of a laboratory connection, the researcher connections become more important 
especially with block 4. In block 4, the links between researchers are strengthened given 
a connection between their laboratories. Researchers in block 3 seem to be the least 
affected by the connections provided by their laboratories. The case of block 1 is quite 
peculiar since it contains two researchers only. This clustering demonstrates that not all 
researchers benefit on equal terms from the institutional level. Some researchers are more 
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P(R=D 



Figure 1: Marginal probabilities of Researcher connections between and within blocks 
(top) and probabilities of Researcher connections between and within blocks condition¬ 
ally on absence (bottom left-hand-side) or presence (bottom right-hand-side) of Lab 
connection. Vertex size is proportional to the block size. Edge width is proportional 
to the probabilities of connection; if this probability is smaller than 0.1, edges are not 
displayed. 


dependent on their laboratories in terms of connections. Furthermore, Figure 2 shows 
the likeliest connections between laboratories are mainly with block 1. The fact that 
researchers are sharing advice inflates the probability of exchanging resources between 
laboratories. 

More importantly, some complex features of the within-level network structure are 
explained mainly by cross-level interactions. In this empirical case, units of each level 
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P(L=1) 


P(L=1|R=0) 




Figure 2: Marginal probabilities of Lab connections between and within blocks (top) 
and probabilities of Lab connections between and within blocks conditionally on absence 
(bottom left-hand-side) or presence (bottom right-hand-side) of Researcher connections. 
Vertex size is proportional to the block size. Edge width is proportional to the proba¬ 
bilities of connection; if this probability is smaller than 0.1, edges are not displayed. 


are clustered in blocks that make sense from the perspective of the categories of people 
and organizations, at each level separately and together. 

Block 1 members work in the biggest labs in terms of size (Figure 3(b)). SBM 
clusters them because they have many more relations than the other members of the 
network. They have among the highest indegrees, and average outdegrees in labs that 
have highest average indegrees and outdegrees (Figures 3(e), 3(f), 3(g), 3(h)). They 
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Table 1: (a) Size of blocks obtained by multiplex SBM. (b) Cross frequencies of blocks 
versus the lab location (in Ile-de-France (IdF) or not), (c) Cross frequencies of blocks 
versus the researcher’s status (lab director or not), (d) Cross frequencies of blocks versus 
the researcher’s specialties (PH: Public Health;Su: Surgery; He: Hematology; ST: Solid 
Tumours; FPh: Fundamental pharmacology; FMo: Fundamental molecular research; 
FMoG: Fundamental molecular genetic research). 



notidf 

idf 

1 

1 

1 

2 

26 

22 

3 

10 

9 

4 

11 

15 



not director 

director 

1 

1 

1 

2 

21 

27 

3 

12 

7 

4 

12 

14 


block 

1 

2 

2 

48 

3 

19 

4 

26 



PH 

Su 

He 

ST 

FPh 

FMo 

FMoG 

1 

0 

0 

0 

0 

0 

1 

1 

2 

12 

5 

6 

10 

7 

7 

1 

3 

0 

1 

3 

1 

1 

10 

3 

4 

6 

1 

7 

4 

2 

2 
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have the highest performance in terms of publication performance scores in both periods 
(Figures 3(c) and 3(d)). In fact, they have a similar relational profile since they are 
providers of transgenic mice for the experiments of many similar colleagues. 

Block 2 members are among the lowest indegrees and outdegrees in labs that have 
the lowest average indegrees and outdegrees (Figures 3(e), 3(f), 3(g), 3(h)), slightly 
older than the others (Figure 3(a)), mainly in the smallest labs in terms of size. This 
block is heterogeneous in terms of specialties (especially 40% clinicians and 25% diag¬ 
nostic/prevention/epidemiology specialists) except fundamental research (Table 1(d)). 
They also have among the lowest performance levels for both periods (Figures 3(c), 3(d)), 
although this is increasing. This is the biggest block. Their behavior may be described 
as fusional as proposed in Lazega et al. (2008) since it corresponds to individuals for 
whom the probabilities of connection are the most affected by the connections of their 
laboratories. 

In block 3, younger fundamental researchers in laboratories carrying out fundamental 
research (mostly molecular research) prevail (Figure 3(a) and Table 1(d)). They have 
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relatively low indegrees and outdegrees, in labs that have the highest indegrees (after 
Block 1) and average outdegrees (Figures 3(e), 3(f), 3(g), 3(h)). 70% are among the top 
performers of this population, i.e. highest performance levels after Block 1 members, for 
both periods (Figures 3(c), 3(d)). Their dominant relational strategies are individualist 
or independentists (as shown in Figure 1) since the probabilities of connection between 
researchers remain quite unchanged, no matter if their laboratories are connected. 

Block 4 is also heterogeneous in terms of specialties but its largest subgroup is com¬ 
posed of hematologists (clinical and fundamentalists) (Table 1(d)). Researchers have 
average indegrees and outdegrees in laboratories that have relatively low indegrees and 
outdegrees and that are also of average size (Figures 3(e), 3(f), 3(g), 3(h)). There are 
proportionally more directors of laboratories in this block than in the others (Table 1(c)). 
Their dominant strategy can be called fusional or collectivist. The performance levels of 
the majority for both periods are somewhat mixed and average, but decreasing (Figures 
3(c), 3(d)). In this block, researchers can take advantage of their laboratories to be 
connected with colleagues but they can also hold relations apart from their laboratories. 

In this case, SBM highlight in the data a specific kind of block structure that provides 
an interesting understanding of the effect of dual positioning. SBM mixes people and 
laboratories that previous analyses used to separate. It is interesting to notice, for 
example, that SBM partitions the population regardless of geographic location (Table 
1 (b)) -a criterion that was previously shown to be meaningful to understand collective 
action in this research milieu. In fact this partition may highlight the links between 
laboratories and between individuals that cut across the boundaries and remoteness 
created by geography, showing that certain categories of actors tend to reach across these 
separations when it suits them, either as investments to prepare future collaborations, 
or as a follow-up for past investments. 


15 





cluster cluster 


(g) (h) 



1 2 3 4 1 2 3 4 

cluster cluster 


16 

Figure 3: Boxplots of age (a), lab size (b), performance in 1st period (c), performance 
in 2nd period (d), indegree for researcher relationship (e), outdegree for researcher re¬ 
lationship (f), indegree for lab relationship (g) and outdegree for lab relationship (h) 
grouped by blocks. 




































4 Discussion 


The essence of ‘networks’ is to help actors cut across organizational boundaries to create 
new relationships (Baker, 1992), to identify new opportunities and, eventually, to create 
new organizations to use or hoard these new opportunities (Lazega, 2012; Lazega et al., 
2013; Tilly, 1998). This work is a step forward a more precise comprehension of such 
mechanisms. 

In this paper, we proposed a SBM modelisation for multiplex networks. This way 
of partitioning the graph and calculating the probabilities that two individuals are con¬ 
nected takes into account more parameters in the graph than just centrality scores and 
size as in Lazega et al. (2007). Besides, compared to MERGMs (Wang et al., 2013), this 
method takes inhomogeneities of actors more into account because the separate blocks 
tend to make sense as blocks with a specific identity. Multiplex SBM statistical analysis 
designed for the analysis of multilevel networks given as “superposed” (Lazega et ah, 
2007) network data aims at identifying the rules that govern the formation of links at 
the inter-individual level conditioned by the characteristics of the inter-organizational 
network. Our analysis found a general structure that shows the ways in which these 
rules can be derived in this case and help assess actor inhomogeneities in terms of their 
influence on parameter estimates. This assessment is new in the sense that it shows that 
formation of connections in the inter-individual network does not apply to all actors in 
the network identically. 

In this paper, the model is applied on K = 2 networks. Applying with K much larger 
would imply estimation difficulties since the number of parameters exponentially depends 
on K. However, particular assumption (such as Markovian dependencies) can reduce the 
difficulty and lead to feasible cases. Besides, in our example, going from multilevel to 
multiplex was quite natural since few laboratories contained more than one researcher. 
If the number of organizations were far smaller then the number of individuals than 
maybe it would be interesting to introduce an other kind of relation, specifying if the 
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individuals are in the same laboratory or not. In this case, the probabilities vector would 
have a special structure, taking this specificity into account (if the individuals belong to 
the same organization then they are automatically linked through their organization) . 

In this version of the analysis, the covariates were studied a posteriori, the classifi¬ 
cation being purely done on the network. From a practical point of view, including the 
covariates in the model would imply estimation difficulty (and even more if the effect of 
the covariates on the connexion probabilities depends on the blocks). However, even if 
attractive at first sight, the choice of including the covariates in the modeling has to be 
questioned. Including them will cluster the data beyond the effect of these covariates, 
while not including them will provide a description of blocks on the basis of the relation, 
may this relation be influenced by these covariates. This second strategy may lead to 
more interpretable results. 

The first aim of Lazega et al. (2008) was cancer research management. Thus, the 
covariate “performance in terms of publication” catches more attention. Indeed, this 
modelisation can generate other questions such as the influence of the networks on the 
performance of the actors (publications in our case). In this work, the performance is 
treated as a factor to explain the relations. However, one could dream of a model which 
would help the actors build the ideal network to optimize the performance. If it is true 
that contemporary society is an “organizational society” (Coleman, 1982; Perrow, 1991; 
Presthus, 1962) -in the sense that action and performance measured at the individual 
level strongly depend on the capacity of the actor to construct and to use organizations 
as instruments, and thus to manage his/her interdependencies at different levels in a 
strategic manner-, then the study of interdependencies jointly at the inter-individual 
and the inter-organizational level is important for numerous sets of problems. Proposing 
a hierarchical model in this direction is out of the scope of the paper but is clearly a 
possible extension. 

To conclude, application to sociological network analysis will clearly benefit from this 


18 


methodological contribution. In spite of several limitations listed before, this analysis 
of multiplex networks seems therefore adapted to certain types of questions that soci¬ 
ologists ask when they try to combine both individual and contextual factors in order 
to estimate the likelihood of an individual or a group to adopt a given behavior or to 
reach a given level of performance. More generally, this approach explores a complex 
meso-social level of accumulation, of appropriation and of sharing of multiple resources. 
This level, still poorly known, is difficult to observe without a structural approach. 
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A Identifiability of the multiplex SBM and convergence of 


the maximum likelihood estimates 
A.l Identifiability 

Celisse et al. (2012) have proved the identifiability of the parameters for the uniplex 
Bernoulli Stochastic Block models. We extend their results to the multiplex SBM. First 
we recall the notations : V(*, j) G {1,. .., n} 2 , i / j, \/(q, l) G {1,..., Q} 2 ,\/w G {0, 1} K , 
we set: 


^ = w\Z i = q,Z J =l) 

a q = P(Zi = q) 

Let us introduce = {^J ) )( q ,i)e{i...Q} and tt = (tt^)^^*. a = (aq,..., a Q ). 
Theorem A.l sets the identifiability of 9 = (a, 7 r) in the multiplex SBM. 

Theorem A.l. Let n > 2 Q. Assume that for any q G {1,... Q}, a q > 0 and for every 
w G {0, 1} K , the coordinates of r^ w \0) = 7 rW'ft are distinct. Then, the multiplex SBM 
parameter 6 = (a,n) is identifiable. 

Proof. The proof given by Celisse et al. (2012) can be directly extended to our case 
leading to the expected result. As explained in the original paper, the identifiability can 
be proved using algebraic tools. The only difference in the multiplex context is that the 
proof has to be applied for any type of tie w G {0,1} A . 

For any w G {0,1}^, we set r q w \o) the probability for a member of block q to have 
a tie of type w with an other individual: r q w \6) = Yl?=i 7T ^ q i' >a i- Let R^ w \6) be the 
Q-square matrix such that R^\o) = {rq‘\o)) 1 for i = 0... Q — 1 and q = 1...Q. 
R (w ) is a Vandermonde matrix, which is invertible by assumptions on the coordinates of 
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r^ ( 9 ). Now, for any w and any i = 0 ... 2 Q — 1, we set: 


Q 

u^\e) = P(X 1 1 i A ' = • ■ • = X]f = W -0) = Y, a M w) (6)Y > ( 6 ) 

q =i 

and M^ w \6) is a (Q + 1) x Q matrix such that 

M^\e) = u < H](6), i = 0...Q,j = 0...Q — 1. (7) 

For any i = 0,. .., Q we define the Q-square matrix M^ w \d) by removing line i from 
this matrix. In particular, 


M QM ^ = R (w) ( 0 ) A ( e ) R M ( 6 y ) ( 8 ) 

where A(6) is the c*-diagonal matrix. All the a q being non-null and R( w \6) being 
invertible, then det> 0. So, if we define: 

Q 

B(X; 6) = ^(-l) i+Q det (M i ^ w \d))X i . 
i =o 

B is of degree Q. Let us define V^ w \d) = (1, ..., (r- w \9))®), then 

B(rj w) (9); 9) = det (lW»(6>), ^ H (0)) , 

where is a (Q + 1) square matrix. The columns of M being linear 

combinations of the V^ w \ we obtain B(r\ w \9); 9) = 0 for any i = 0... Q — 1. So, we 
can factorize B as 

Q -1 

B(x; 9) = det (M Q ^ {9)) JJ (x - r\ w) ( 9 )). (9) 

i =0 

Now assume that 9 = ( 7 r,a) and = (7r / ,o: / ) are two sets or parameters such that 
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for any multiplex graph X 1;A , £(X 1;A ;$) = £(X. 1 :K -,9'). Consequently, from equation 
(6), we get u[ w \d) = u\ w \d') and from equation (7), we deduce M l ( w \ 8 ) = MA W )(Q') 
for any i = 0,... ,Q — 1. The polynomial function B depending on the determinant 
of the M^ w \9)’s, we have B(-;9 ) = B(-\9'), leading to ri(w){9) = ri(w)(9’) for all 
i = 0,..., Q — 1 thanks to equation (9). Thus R^ w \9) = R^ w \9') and 

A(9) = (R^\9y)- 1 M Q ^ w \9)R iw '>{9) = A(9 '). 

As a consequence, a = a'. 

Finally, let (0 > i,j , > Q) denote : 

= P (Xl'j? = w, k = 1 ,..., i + 1 and k = n — j + 1,..., n). 

We can write the Q x Q matrix U^ w \9) as C/H(0) = R^ w \9)A(9)7v^A{9){R^ w \9)y. 
Using the fact that R^(9) = R^ w \ 6 ') and A{9) = A{ 6 '), we get 

UW( 6 ) = U^ w \ 6 ') => = (tt'W), Vw e { 0 , \} K . 

And the theorem is demonstrated. 

□ 

A.2 Consistency of the maximum likelihood estimator in multiplex 

SBM models 

We study the asymptotic properties of the maximum likelihood estimator in multiplex 
SBM models. Let Ai, A 2 and A 3 be the following assumptions: 

Assumption Ai : for every q A q', there exists l E {1,..., Q} such that n q 1 / 7 v q q , or 
T^lq Itlq' • 

Assumption there exists C, > 0 such that \/(q, l ) E 1,..., Q, 1 E]0, l[=i> ir^ E 
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Assumption A 3 : there exists 7 E (0,1 /Q) such that Vg G 1, ..., Q, a q G [7, 1 — 7]. 


Let X 1:A be a realization from the multiplex SBM: P( • ] Z*, a*, iv*) where Z* = 
(Z (,..., Z*) is the true group label sequence and (a*, 7r*) are the true parameters. Let 
(S,7r) be the maximum likelihood estimator defined as: 


(S, 7 r) = Argmax (a 7r) £(X 1:Ar ; a, tt) , 


where 

n 

e 

Ze{i,...Q} n t=l 

n n 

E logP(X^|Z i ,Z j ;7r)= E E ^WogTrJE 

wAj iue{0,l} x 

Theorem A.2. Let (A\), (A 2 ) and (A 3 ) hold. Then, for any distance d(-,-) on the set 
of parameter tv , we have: 

d(n, tv*) —-— > 0 . 
n—>00 

Moreover, assume that ||7T — 7T* Hoc = op(y / log(n) / n) then, for any distance d(-, •) in 

d(a, a*) —— > 0 . 

n—> 00 

Note the rate op(-y/log(n)/ n) has not been proved yet. However, in the unilevel 
context, there is empirical evidence that the rate convergence on fv is 1/n (Gazal et ah, 
2012 ). 


i(X 1 :K -a, tv) = 
£i(X\Z;tv) = 
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B Variational EM algorithm for multiplex SBM : princi¬ 
ple, details and convergence 

B.l General principal of the variational EM 

The Stochastic Block Models belong to the incomplete data models class, the non¬ 
observed data being the block indices (Zi)i= i... n 6 {1,. .. ,Q} n - As written before, the 
likelihood has a marginal expression: 

£(X 1 :K -,9) = e ^ i ( xl:iC l Z; 7 r )p(z ; a ), (10) 

ze{i ,...Q} n 

which is not tractable as soon as n and Q are large. The variational EM is an alternative 
method to maximize the marginal likelihood with respect to 6 . The variational EM (ap¬ 
plied in the SBM context by Daudin et al. (2008)) relies on the following decomposition 
of (10). Let 1Z xi -K be any probability distribution on Z, wg havoi 

log^(X 1:A >) = ^^ xl; K(Z)logp(X 1 :K ,Z;0)-^^ xl :K(Z)log^ x i : K(Z) 

z z 

+KL[ft x i=*,p(-|X 1 :A: ;0)], ( 11 ) 

where KL is the Kullback-Leibler distance, p(X l:K ,Z;9) is the joint density of X 1:A 
and Z (namely the complete likelihood) and p(-\X 1:K ;9) is the posterior distribution of 
Z given the data X 1:A and the parameters 6 . Instead of maximizing logt?(X 1:A ; 9), the 
variational E1V1 optimizes a lower bound Hq (7Z X i, k ) of log^(X 1:A |0) where: 

MK* i:k) = log{(X 1:K l9)-KL[7l x i:K,p(-lX hK ;9)], 

= Y 7Z x^(Z)logp(X 1:K ,Z;9)~Y 7Z x^(Z)log7Z xl :K(Z), (12) 

z z 

< log£(X 1 :K l9). 
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Note that, thanks to equality (12), optimizing Xq{1Z-^i-.k) with respect to 6 no longer re¬ 
quires the computation of the marginal likelihood. Note also that the equality Xq(JZ-^i-.k) = 
log £(X. 1:K -,0) holds if and only if 1Z -K — p( • | ^ 5 @) • As ci consequence , .k will 

be taken as an approximation of p(-|X 1: ^; 6) in a certain class of distributions. Jaakkola 
(2000) proposed to optimize it in the following class: 

n 

= P h(Zi,Ti ), 
i =1 

where /i(-;Tj) is the multinomial distribution of parameter = (r*i,... ,Tj g ). Finally, 
the variational EM updates alternatively 6 and r in the following way. At iteration (f), 
given the current state 

- Step 1 Compute 

= argmin T KL[7^x 1: Jf iT ,p(-|X 1:ft ‘; = arg max T X^t-p (7il x i:K T ). 

- Step 2 Compute 6^ = argmax^X^p T (t))- 

The details of steps 1 and 2 directly depend on the considered statistical model. For 
uniplex SBM without covariates, they are given in Daudin et al. (2008). The details for 
the multiplex SBM are given here after. 

B.2 Details of the calculus for multiplex SBM models 

We now detail Step 1 and Step 2 for multiplex SBM models. 

• Step 1: T 0) verifies 

= argminKL[7?, x i ; x T ,p(-\X. 1:K -, 0^ -1 ))] = argmaxL(i-i)(K x i:K T ). 

T ’ T ’ 
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We first rewrite Iq(TZ x i-.k t ) for this special context: 


le{TZ x i.K T ) 


Y ^x«,r( Z ) logp(X 1:Ar , Z; 0) - £ ft x i=* T (Z) log 7^ x i=x T (Z), 

z z 


with 


l°gp(X 1:K , Z; 0) = f?i(X i;A |Z; 0) + logp(Z; 9 ), 

n 

= Y log P( X ij K IX*, Zj ; 0) + Y !°g ■ 

hiMi i=1 

This quantity has to be integrated over Z where Z ~ 1Z x v.k t which means that Z 
(Zi)i= i... n are independent variables such that P(Zj = q) = Ti q . We obtain: 


1b(K x i : k >T ) = Y' R '^. 1 '- k ,t{ Z 

z 

-Y^"«A Z )^K X1 :K, t (Z) , 


X] lo g \Zi,Zj-9) + Y lo & a Zi 


2=1 


Y l °SP( X ij K \Zi = 9, Zj = l; e)T iq Tji 
q,i 

n Q n Q 

+EE T iq log a q -EE Tig log T iq , 

i=l g=l i= 1 g=l 


where logp(ATb A |Zj = q, Zj = l; 0)’s expression is given in equation (3). 

Zo(TZ x i-,K T ) has to be maximized with respect to t under the constraint: Vi = 
1 • • • n > E?=i — 1. A.S cl CORSGQUGILCG, WG COHipiltG tllG d.Griva,tivGS of X (9 T ) + 

Eti A. E?.i % - 1 with respect to (Ai)i=i... n and {Ti q )i = i... ntq =i...Q where Ai are the 
Lagrange multipliers, leading to the following collection of equations: for i = 1... n and 
q = l...Q, 
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n 

^2 ^2 l °&P( X ij K \ Zi = < h Z j = l '> °) T jl + lo S a q ~ lo S T iq + l + Xi=0, 
l j=l> j¥=i 

which leads to the following fixed point problem: 


n Q 

Ti q = e 1+Xi a q Y\p(X}j K \Zi = q,Zj = l;9) T i l , Vi = 1... n,\/q = 1... Q , 

j= IjV* t —i 


which has to be solved under the constraints Vi = 1... n, T iq = T This optimiza¬ 

tion problem is solved using a standard fixed point algorithm. 

Step 2 Compute 9^ = arg max# T e ( t ) (TZ-^uk T (t))- 

Once the f have been optimized, the parameters 8 maximizing Ig(TZ^i-.K f) have 
to be computed under the constraints: Y2^=i a q = 1 and Ylwe (o,i} l = 1 for a11 

(q, 0 £{!,•• • ,Q} 2 - 

The maximization with respect to a is quite direct and in any case, we obtain: 


1 

a i = 

1=1 


TT iw) ~ 

n qi ~ 


yii.j TiqTj]I X l:K =ru 

Thii T iq T ji 


Remark 5. If the edge probabilities depend on covariates: 


logit{TT ( q 2) = 2 q 2 + (P^Yyij , 

then the optimization of (2^) am ^ (P^) a t s ^ e P ® of the VarEm is not explicit anymore 
and one should resort to optimization algorithms such as Newton-Raphson algorithm. 
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B.3 Convergence of the VarEM estimates 


We now consider the consistency of the estimates obtained by the variational EM algo¬ 
rithm. Using the variational EM previously described is equivalent to maximizing the 
so-called variational-likelihood Z where: 

X(X i:A ; T, a, 7T) = T iqTjl ^ hx ij= w log TV^ - ^ Tj g (log T iq ~ log Oi q ) . 

bU w (o,i} K ig 

The variational estimators (VE) are obtained by: 

r(o:,7r) = argrnax X(X 1:A ‘; t, a, 7r) (5, 7 t) = argmax Z(X 1:A q t, a, 7 t) 

T CX..7T 

Theorem B.l. Assume that (A\), (A 2 ) and, (A 3 ) hold. Then for any distance on the 
set of parameters tv, 

d(TV, 7T*) —)> 0 . 

n—>00 

Moreover, assume that d( tv, tv*) = op(l/ro), then for any distance on 

d(a, a*) —0 . 

n—>■ 00 

B.4 About the proofs 

The proofs of these results require many intermediate results which won’t be given here 
because their adaptation to the multiplex context is quite direct. Indeed, at any step 
of the proof, the original Bernoulli distribution needs to replaced by its A'-dimensional 
version. More precisely, a central quantity in the proof is the ratio log . In 
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the unilevel case, this quantity is: 


log 


P(Z|X i:A ; a, 7r) 
p( Z*|X 1:A ;a,7r) 


= lo8 flX^|Z;.) p(Z;a) 


>(X 1:A |Z*;7r) 
(w) 


' p(Z*; a) 


7 r 


X x « lo s-£r + <1 “ x « )log 

*A? ^ 


1 — 7T 


(jd) 
Z-i, Zj 


7T ( ™K 

z . ,z* 

1 ’ J 


i - VI 


+ 


n 

x 

i= 1 


i a Zi 

log—i 
a z * 


In the multiplex case, the sum of two terms is replaced by a sum of 2 K terms: 


log 


p(Z|X;a,7r) 

p(Z*|X; a, 7r) 


= X X 

*A? «>e{o,i} K 


n ( w ) 

Vo-} ^ + 


1 a Zi 

log- 


7T 


z*,z* 


i= 1 


CK, 


where S«;e{o i} K 7 T qi’^ = 1> f° r an y (<b0- Going from two terms to 2 A terms does not 
imply any mathematical difficulty and so does not compromise the convergence results. 
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